OHI-Northeast | OHI Science | Citation policy
For OHI Global Assessments, we calculated a climatology from 1982-2017 per cell using weekly data from CorTAD. We then calculated the number of times a given cell’s weekly SST was greater than the climatological mean for that week (an anomaly: greater than mean + one standard deviation) and summed the number of weekly anomalies in a single year. The maximum value a cell could have is 52 which would mean that cell had anomalous SST temperatures for each week of the year.
To account for annual variation, we look at Sea Surface Temperature anomalies in 5 year periods, so the maximum value possible per cell is 260 anomalous weeks. To rescale the values from 0 to 1 we set a reference point. Previously, the reference point for SST has just been the maximum difference in anomalous weeks between the most recent time period and a historical reference period (1985-1989).
This time we have decided to use a reference point that represents a regime shift. Once a given cell is anomalous for more than 50% of a five-year period, it has shifted into a new regime. All cells that have a value greater than 130 weeks (51% of a 5 year time period) are assigned a 1. The rest of the cells are scaled to this reference point by dividing by 130.
Source: The Coral Reef Temperature Anomaly Database Version 6 (CoRTAD)
Downloaded: August 21, 2018
Description: Sea Surface Temperature Anomalies (Kelvin)
Native data resolution: 4km2
Time range: 1982 - 2017
Format: NetCDF
We’re going to use the global data that was processed for OHI 2018. This data is held on a server at NCEAS.
Each of these rasters is the number of positive weekly anomalies each year.
sst_global_files <- list.files(file.path(dir_M,'git-annex/globalprep/prs_sst/v2018/int'), pattern = 'annual_pos_anomalies', full.names=T)
plot(raster(sst_global_files[36]),col=cols,main = "Sea Surface Temperature 2017",box=F,axes=F,
legend.args=list(text='Anomalous Weeks', side=4, font=2, line=2.5, cex=0.8))Using the crop function from the raster package we crop all sea surface temperature rasters to our extent and then reproject them to the US Albers projection for consistency across the assessment. We crop the global rasters first to reduce the time it takes to reproject the rasters. ocean_ne is used as a mask to remove land cells from the raster for better visuals.
registerDoParallel(10)
foreach(f = sst_global_files) %dopar% {
raster(f)%>% #raster the file
crop(wgs_ext)%>% #crop to the WGS extent
projectRaster(ocean_ne)%>%
mask(ocean_ne, filename = paste0(file.path(dir_anx),'/prs_sst/output/sst_annual_anoms/annual_anoms_', substr(basename(f), 25, 28),'.tif'), overwrite = T)
}
plot(raster(file.path(dir_anx, 'prs_sst/output/sst_annual_anoms/annual_anoms_2017.tif')),col=cols, axes=F, main = "Sea Surface Temperature 2017",
legend.args=list(text='Anomalous Weeks', side=4, font=2, line=2.5, cex=0.8))Getting a sense of how things have changed over time.
#rasterize each file in the list of files and then create a brick of rasters
s <- lapply(list.files(file.path(dir_anx, 'prs_sst/output/sst_annual_anoms'), full.names=T), raster) %>%
brick()
names(s) <- paste0("Years_",(substr(names(s), 14, 17))) #rename each layer for plotting
gsub(".", "-", names(s), fixed = TRUE) #replace the . with a -
library(animation)
saveGIF({
for(i in 1:nlayers(s)){
# don't forget to fix the zlimits
plot(s[[i]], zlim=c(0,52), axes=F, col=cols,
main=names(s[[i]]),
legend.args=list(text='Anomalous Weeks', side=4, font=2, line=2.5, cex=0.8))
}
}, movie.name = 'sst_annual_anoms.gif')Calculating total anomalous weeks for each 5-year period from 1982 - 2017.
l <- list.files(file.path(dir_anx, 'prs_sst/output/sst_annual_anoms'), full.names=T)
for(i in 1982:2013){ #i=2005
yrs <- c(i, i+1,i+2, i+3, i+4)
s <- stack(l[substr(l, 80, 83) %in% yrs]) %>%
sum(.)
writeRaster(s, filename = paste0(dir_anx, '/prs_sst/output/sst_5yr_anom_sums/sum_anoms_', min(yrs), '-', max(yrs), '.tif'), overwrite=T)
}To account for annual variation, we look at Sea Surface Temperature anomalies in 5 year periods, so the maximum value possible per cell is 260 anomalous weeks. To rescale the values from 0 to 1 we need to set a reference point. Previously, the reference point for SST has just been the maximum difference in anomalous weeks between the most recent time period and a historical reference period (1985-1989).
This time we have decided to use a reference point that represents a regime shift. Once a given cell is anomalous for more than 50% of a five-year period, it has shifted into a new regime. All cells that have a value greater than 130 weeks (51% of a 5 year time period) are assigned a 1. The rest of the cells are scaled to this reference point by dividing by 130.
sst_aggs <- list.files(file.path(dir_anx, 'prs_sst/output/sst_5yr_anom_sums'), full.names=T)
resc_func <- function(x){
#get the year from the file for naming the output raster
yrs <- substr(x, 78, 86)
#if a cell value is greater than or equal to the reference point (130 weeks), we assign a value of 1, otherwise it is divided by the reference point
raster(x)%>%
calc(., fun=function(x){ifelse(x < 0, 0, ifelse(x > 130, 1, x/130))},
filename = paste0(dir_anx, '/prs_sst/output/sst_rescale/sst_rescale_', yrs, '.tif'), overwrite=T)
}
foreach(file = sst_aggs) %dopar%{
resc_func(file)
}Using the local reference point, we can see how the pressure changes over time.
#list the rescaled rasters and assign them in a brick
resc <- lapply(list.files(file.path(dir_anx, 'prs_sst/output/sst_rescale'), full.names=T), raster) %>%
brick()
names(resc) <- paste0("Years_",(substr(names(resc),13,21)))
gsub(".", "-", names(resc), fixed = TRUE)
#create a gif of values over time
saveGIF({
for(i in 1:nlayers(resc)){
# don't forget to fix the zlimits
plot(resc[[i]], zlim=c(0,1), axes=F, col=cols,
main=names(resc[[i]]),
legend.args=list(text='Anomalous Weeks', side=4, font=2, line=2.5, cex=0.8))
}
}, movie.name = 'sst_rescale.gif')By extracting the data for each of the 9 regions using the zonal function from the raster package we can get the mean score per region.
# read in rescaled files
pressure_stack <- lapply(list.files(file.path(dir_anx, 'prs_sst/output/sst_rescale'),full.names=T), raster)%>%
brick()
# extract data for each region (zones is previously defined in common.R)
regions_stats <- zonal(pressure_stack, zones, fun="mean", na.rm=TRUE, progress="text")%>%
data.frame()
#merge the zonal stats with the rgn_data as defined in common.R
data <- merge(rgn_data, regions_stats, all.y=TRUE, by.x="rgn_id", by.y="zone") %>%
dplyr::select(-area_km2)%>%
gather("year", "pressure_score", starts_with("sst_rescale_"))
#add a year column that grabs the last year of the 5 year period and reformat as numeric
sst_data <- data %>%
mutate(year=substr(year,18, 21)) %>%
mutate(year = as.numeric(year))
#save output data
write.csv(sst_data, file.path(dir_calc, "layers/cc_sst.csv"), row.names=FALSE)Let’s look at the most recent scores using the data from 2012 since it’s the most recent year in the CoRTAD dataset.
#select data from 2017
now <- read.csv(file.path(dir_calc, "layers/cc_sst.csv")) %>%
filter(year==2017)
#map_scores is a function to plot a tmap map of the scores
map_scores(now, score_var = now$pressure_score,
scale_label = "Pressure Score",
map_title = "Sea Surface Temperature Pressures 2017")library(ggplot2)
sst_data <- read_csv(file.path(dir_calc, "layers/cc_sst.csv"))
ggplot(sst_data, aes(x = year, y = pressure_score, color = rgn_name))+
geom_line()+
labs(color = "Region") +
theme_bw()This code chunk creates a google visualization of the scores through time. It unfortunately does not allow the google chart to show up in the RMarkdown knitted document.
library(googleVis)
plotData <- sst_data %>%
dplyr::select(rgn_name, year, pressure_score)
Motion=gvisMotionChart(plotData,
idvar="rgn_name",
timevar="year",
options=list(width=550, height=450))
print(Motion, 'chart', file='sst.html')Saha, Korak; Zhao, Xuepeng; Zhang, Huai-min; Casey, Kenneth S.; Zhang, Dexin; Zhang, Yongsheng; Baker-Yeboah, Sheekela; Relph, John M.; Krishnan, Ajay; Ryan, Thomas (2018). The Coral Reef Temperature Anomaly Database (CoRTAD) Version 6 - Global, 4 km Sea Surface Temperature and Related Thermal Stress Metrics for 1982 to 2017. NOAA National Centers for Environmental Information. Dataset. https://doi.org/10.25921/ffw7-cs39. Accessed: August 21, 2018/
Selig, E.R., K.S. Casey, and J.F. Bruno (2010), New insights into global patterns of ocean temperature anomalies: implications for coral reef health and management, Global Ecology and Biogeography, DOI: 10.1111/j.1466-8238.2009.00522.x.